Terrain Analysis for Drone Survey Areas

The original DTM for the survey area was produced in Pix4D from the RGB image data collected by the DJI P4 Multispectral Drone Surveys. Terrain analysis was carried out using the terrain function from the terra package.

Bokspits 1 Survey Area

##----1. Terrain analysis Bokspits 1 survey area----
##Import shape file data

Bokspits_1_clipper <- read_sf(dsn = 'E:/Glenn/Botswana/Final_Drone_Survey_Data/Bokspits_1', layer = "Bokspits_1_clip")
Bokspits_1_clip <- vect(Bokspits_1_clipper)

##Import raster DTM File
Bokspits_1_DTM <-  rast("E:/Glenn/Botswana/Pix4d/Bokspits_1_MS_RGB/3_dsm_ortho/extras/dtm//Bokspits_1_MS_RGB_dtm.tif")

plot(Bokspits_1_DTM)

Slope in degrees

##terrain(x, opt='slope', unit='radians', neighbors=8, filename='', ...) - NB This is from the original documentation but it doesnt work - need to omit opt=
Bokspits_1_DTM_slope <- terrain(Bokspits_1_DTM,'slope', unit='degrees')
plot(Bokspits_1_DTM_slope)

Aspect in degrees

Bokspits_1_DTM_aspect <- terrain(Bokspits_1_DTM,'aspect', unit='degrees')
plot(Bokspits_1_DTM_aspect)

Surface roughness

Bokspits_1_DTM_roughness <- terrain(Bokspits_1_DTM,'roughness')
plot(Bokspits_1_DTM_roughness)

Surface water flow direction

Bokspits_1_DTM_flowdir <- terrain(Bokspits_1_DTM,'flowdir')
plot(Bokspits_1_DTM_flowdir)

Individual data sets were combined into a multilayer raster stack of terrain features and cropped to the same extent as the calibrated reflection image data for Bokspits 1.

## Make image stack of terrain analysis layers cropped to study area boundary
Bokspits_1_Stack_Terrain <- c(Bokspits_1_DTM_slope,Bokspits_1_DTM_aspect,Bokspits_1_DTM_roughness,Bokspits_1_DTM_flowdir ) 
Bokspits_1_Stack_TerrainCr <- crop(Bokspits_1_Stack_Terrain,Bokspits_1_clip )
Bokspits_1_Stack_TerrainCrop <- mask(Bokspits_1_Stack_TerrainCr,Bokspits_1_clip )
plot(Bokspits_1_Stack_TerrainCrop)

Bokspits 2 Survey Area

Bokspits_2_clipper <- read_sf(dsn = 'E:/Glenn/Botswana/Final_Drone_Survey_Data/Bokspits_2', layer = "Bokspits_2_clip")
Bokspits_2_clip <- vect(Bokspits_2_clipper)

##Import raster DTM File
Bokspits_2_DTM <-  rast("E:/Glenn/Botswana/Pix4d/Bokspits_2_MS_RGB/3_dsm_ortho/extras/dtm//Bokspits_2_MS_RGB_dtm.tif")

plot(Bokspits_2_DTM)

Slope

Bokspits_2_DTM_slope <- terrain(Bokspits_2_DTM,'slope', unit='degrees')
plot(Bokspits_2_DTM_slope)

Aspect

Bokspits_2_DTM_aspect <- terrain(Bokspits_2_DTM,'aspect', unit='degrees')
plot(Bokspits_2_DTM_aspect)

Roughness

Bokspits_2_DTM_roughness <- terrain(Bokspits_2_DTM,'roughness')
plot(Bokspits_2_DTM_roughness)

Surface water flow direction

Bokspits_2_DTM_flowdir <- terrain(Bokspits_2_DTM,'flowdir')
plot(Bokspits_2_DTM_flowdir)

Individual data sets were combined into a multilayer raster stack of terrain features and cropped to the same extent as the calibrated reflection image data for Bokspits 2.

## Make image stack of terrain analysis layers cropped to study area boundary
Bokspits_2_Stack_Terrain <- c(Bokspits_2_DTM_slope,Bokspits_2_DTM_aspect,Bokspits_2_DTM_roughness,Bokspits_2_DTM_flowdir ) 
Bokspits_2_Stack_TerrainCr <- crop(Bokspits_2_Stack_Terrain,Bokspits_2_clip )
Bokspits_2_Stack_TerrainCrop <- mask(Bokspits_2_Stack_TerrainCr,Bokspits_2_clip )
plot(Bokspits_2_Stack_TerrainCrop)

Bokspits 3 Survey Area

Bokspits_3_clipper <- read_sf(dsn = 'E:/Glenn/Botswana/Final_Drone_Survey_Data/Bokspits_3', layer = "Bokspits_3_clip")
Bokspits_3_clip <- vect(Bokspits_3_clipper)

##Import raster DTM File
Bokspits_3_DTM <-  rast("E:/Glenn/Botswana/Pix4d/Bokspits_3_MS_RGB/3_dsm_ortho/extras/dtm//Bokspits_3_MS_RGB_dtm.tif")

plot(Bokspits_3_DTM)

Slope

Bokspits_3_DTM_slope <- terrain(Bokspits_3_DTM,'slope', unit='degrees')
plot(Bokspits_3_DTM_slope)

Aspect

Bokspits_3_DTM_aspect <- terrain(Bokspits_3_DTM,'aspect', unit='degrees')
plot(Bokspits_3_DTM_aspect)

Roughness

Bokspits_3_DTM_roughness <- terrain(Bokspits_3_DTM,'roughness')
plot(Bokspits_3_DTM_roughness)

Surface water flow direction

Bokspits_3_DTM_flowdir <- terrain(Bokspits_3_DTM,'flowdir')
plot(Bokspits_3_DTM_flowdir)

Individual data sets were combined into a multilayer raster stack of terrain features and cropped to the same extent as the calibrated reflection image data for Bokspits 3.

## Make image stack of terrain analysis layers cropped to study area boundary
Bokspits_3_Stack_Terrain <- c(Bokspits_3_DTM_slope,Bokspits_3_DTM_aspect,Bokspits_3_DTM_roughness,Bokspits_3_DTM_flowdir ) 
Bokspits_3_Stack_TerrainCr <- crop(Bokspits_3_Stack_Terrain,Bokspits_3_clip )
Bokspits_3_Stack_TerrainCrop <- mask(Bokspits_3_Stack_TerrainCr,Bokspits_3_clip )
plot(Bokspits_3_Stack_TerrainCrop)

Struizendam 1 Survey Area

##----1. Terrain analysis Struizendam 1 survey area----
##Import shape file data

Struizendam_1_clipper <- read_sf(dsn = 'E:/Glenn/Botswana/Final_Drone_Survey_Data/Struizendam_1', layer = "Struizendam_1_clip")
Struizendam_1_clip <- vect(Struizendam_1_clipper)

##Import raster DTM File
Struizendam_1_DTM <-  rast("E:/Glenn/Botswana/Pix4d/Struizendam_1_MS_RGB/3_dsm_ortho/extras/dtm//Struizendam_1_MS_RGB_dtm.tif")

plot(Struizendam_1_DTM)

Slope in degrees

##terrain(x, opt='slope', unit='radians', neighbors=8, filename='', ...) - NB This is from the original documentation but it doesnt work - need to omit opt=
Struizendam_1_DTM_slope <- terrain(Struizendam_1_DTM,'slope', unit='degrees')
plot(Struizendam_1_DTM_slope)

Aspect in degrees

Struizendam_1_DTM_aspect <- terrain(Struizendam_1_DTM,'aspect', unit='degrees')
plot(Struizendam_1_DTM_aspect)

Surface roughness

Struizendam_1_DTM_roughness <- terrain(Struizendam_1_DTM,'roughness')
plot(Struizendam_1_DTM_roughness)

Surface water flow direction

Struizendam_1_DTM_flowdir <- terrain(Struizendam_1_DTM,'flowdir')
plot(Struizendam_1_DTM_flowdir)

Individual data sets were combined into a multilayer raster stack of terrain features and cropped to the same extent as the calibrated reflection image data for Struizendam 1.

## Make image stack of terrain analysis layers cropped to study area boundary
Struizendam_1_Stack_Terrain <- c(Struizendam_1_DTM_slope,Struizendam_1_DTM_aspect,Struizendam_1_DTM_roughness,Struizendam_1_DTM_flowdir ) 
Struizendam_1_Stack_TerrainCr <- crop(Struizendam_1_Stack_Terrain,Struizendam_1_clip )
Struizendam_1_Stack_TerrainCrop <- mask(Struizendam_1_Stack_TerrainCr,Struizendam_1_clip )
plot(Struizendam_1_Stack_TerrainCrop)

Struizendam 2 Survey Area

Struizendam_2_clipper <- read_sf(dsn = 'E:/Glenn/Botswana/Final_Drone_Survey_Data/Struizendam_2', layer = "Struizendam_2_clip")
Struizendam_2_clip <- vect(Struizendam_2_clipper)

##Import raster DTM File
Struizendam_2_DTM <-  rast("E:/Glenn/Botswana/Pix4d/Struizendam_2_MS_RGB/3_dsm_ortho/extras/dtm//Struizendam_2_MS_RGB_dtm.tif")

plot(Struizendam_2_DTM)

Slope

Struizendam_2_DTM_slope <- terrain(Struizendam_2_DTM,'slope', unit='degrees')
plot(Struizendam_2_DTM_slope)

Aspect

Struizendam_2_DTM_aspect <- terrain(Struizendam_2_DTM,'aspect', unit='degrees')
plot(Struizendam_2_DTM_aspect)

Roughness

Struizendam_2_DTM_roughness <- terrain(Struizendam_2_DTM,'roughness')
plot(Struizendam_2_DTM_roughness)

Surface water flow direction

Struizendam_2_DTM_flowdir <- terrain(Struizendam_2_DTM,'flowdir')
plot(Struizendam_2_DTM_flowdir)

Individual data sets were combined into a multilayer raster stack of terrain features and cropped to the same extent as the calibrated reflection image data for Struizendam 2.

## Make image stack of terrain analysis layers cropped to study area boundary
Struizendam_2_Stack_Terrain <- c(Struizendam_2_DTM_slope,Struizendam_2_DTM_aspect,Struizendam_2_DTM_roughness,Struizendam_2_DTM_flowdir ) 
Struizendam_2_Stack_TerrainCr <- crop(Struizendam_2_Stack_Terrain,Struizendam_2_clip )
Struizendam_2_Stack_TerrainCrop <- mask(Struizendam_2_Stack_TerrainCr,Struizendam_2_clip )
plot(Struizendam_2_Stack_TerrainCrop)

Struizendam 3 Survey Area

Struizendam_3_clipper <- read_sf(dsn = 'E:/Glenn/Botswana/Final_Drone_Survey_Data/Struizendam_3', layer = "Struizendam_3_clip")
Struizendam_3_clip <- vect(Struizendam_3_clipper)

##Import raster DTM File
Struizendam_3_DTM <-  rast("E:/Glenn/Botswana/Pix4d/Struizendam_3_MS_RGB/3_dsm_ortho/extras/dtm//Struizendam_3_MS_RGB_dtm.tif")

plot(Struizendam_3_DTM)

Slope

Struizendam_3_DTM_slope <- terrain(Struizendam_3_DTM,'slope', unit='degrees')
plot(Struizendam_3_DTM_slope)

Aspect

Struizendam_3_DTM_aspect <- terrain(Struizendam_3_DTM,'aspect', unit='degrees')
plot(Struizendam_3_DTM_aspect)

Roughness

Struizendam_3_DTM_roughness <- terrain(Struizendam_3_DTM,'roughness')
plot(Struizendam_3_DTM_roughness)

Surface water flow direction

Struizendam_3_DTM_flowdir <- terrain(Struizendam_3_DTM,'flowdir')
plot(Struizendam_3_DTM_flowdir)

Individual data sets were combined into a multilayer raster stack of terrain features and cropped to the same extent as the calibrated reflection image data for Struizendam 3.

## Make image stack of terrain analysis layers cropped to study area boundary
Struizendam_3_Stack_Terrain <- c(Struizendam_3_DTM_slope,Struizendam_3_DTM_aspect,Struizendam_3_DTM_roughness,Struizendam_3_DTM_flowdir ) 
Struizendam_3_Stack_TerrainCr <- crop(Struizendam_3_Stack_Terrain,Struizendam_3_clip )
Struizendam_3_Stack_TerrainCrop <- mask(Struizendam_3_Stack_TerrainCr,Struizendam_3_clip )
plot(Struizendam_3_Stack_TerrainCrop)

Struizendam 4 Survey Area

##----4. Terrain analysis Struizendam 4 survey area----
##Import shape file data

Struizendam_4_clipper <- read_sf(dsn = 'E:/Glenn/Botswana/Final_Drone_Survey_Data/Struizendam_4', layer = "Struizendam_4_clip")
Struizendam_4_clip <- vect(Struizendam_4_clipper)

##Import raster DTM File
Struizendam_4_DTM <-  rast("E:/Glenn/Botswana/Pix4d/Struizendam_4_MS_RGB/3_dsm_ortho/extras/dtm//Struizendam_4_MS_RGB_dtm.tif")

plot(Struizendam_4_DTM)

Slope in degrees

##terrain(x, opt='slope', unit='radians', neighbors=8, filename='', ...) - NB This is from the original documentation but it doesnt work - need to omit opt=
Struizendam_4_DTM_slope <- terrain(Struizendam_4_DTM,'slope', unit='degrees')
plot(Struizendam_4_DTM_slope)

Aspect in degrees

Struizendam_4_DTM_aspect <- terrain(Struizendam_4_DTM,'aspect', unit='degrees')
plot(Struizendam_4_DTM_aspect)

Surface roughness

Struizendam_4_DTM_roughness <- terrain(Struizendam_4_DTM,'roughness')
plot(Struizendam_4_DTM_roughness)

Surface water flow direction

Struizendam_4_DTM_flowdir <- terrain(Struizendam_4_DTM,'flowdir')
plot(Struizendam_4_DTM_flowdir)

Individual data sets were combined into a multilayer raster stack of terrain features and cropped to the same extent as the calibrated reflection image data for Struizendam 4.

## Make image stack of terrain analysis layers cropped to study area boundary
Struizendam_4_Stack_Terrain <- c(Struizendam_4_DTM_slope,Struizendam_4_DTM_aspect,Struizendam_4_DTM_roughness,Struizendam_4_DTM_flowdir ) 
Struizendam_4_Stack_TerrainCr <- crop(Struizendam_4_Stack_Terrain,Struizendam_4_clip )
Struizendam_4_Stack_TerrainCrop <- mask(Struizendam_4_Stack_TerrainCr,Struizendam_4_clip )
plot(Struizendam_4_Stack_TerrainCrop)

sessionInfo()
## R version 4.1.2 (2021-11-01)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United Kingdom.1252 
## [2] LC_CTYPE=English_United Kingdom.1252   
## [3] LC_MONETARY=English_United Kingdom.1252
## [4] LC_NUMERIC=C                           
## [5] LC_TIME=English_United Kingdom.1252    
## 
## attached base packages:
## [1] splines   stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] rgeos_0.5-9        sp_1.4-7           terra_1.5-21       cowplot_1.1.1     
##  [5] writexl_1.4.0      sf_1.0-7           DescTools_0.99.45  gridExtra_2.3     
##  [9] MASS_7.3-54        RColorBrewer_1.1-3 lubridate_1.8.0    forcats_0.5.1     
## [13] stringr_1.4.0      dplyr_1.0.9        purrr_0.3.4        readr_2.1.2       
## [17] tidyr_1.2.0        tibble_3.1.7       ggplot2_3.3.6      tidyverse_1.3.1   
## [21] viridis_0.6.2      viridisLite_0.4.0 
## 
## loaded via a namespace (and not attached):
##  [1] fs_1.5.2           httr_1.4.3         tools_4.1.2        backports_1.4.1   
##  [5] bslib_0.3.1        utf8_1.2.2         R6_2.5.1           KernSmooth_2.23-20
##  [9] DBI_1.1.2          colorspace_2.0-3   withr_2.5.0        tidyselect_1.1.2  
## [13] Exact_3.1          compiler_4.1.2     cli_3.3.0          rvest_1.0.2       
## [17] expm_0.999-6       xml2_1.3.3         sass_0.4.0         scales_1.2.0      
## [21] classInt_0.4-3     mvtnorm_1.1-3      proxy_0.4-26       digest_0.6.29     
## [25] rmarkdown_2.14     pkgconfig_2.0.3    htmltools_0.5.2    highr_0.9         
## [29] dbplyr_2.1.1       fastmap_1.1.0      rlang_1.0.2        readxl_1.4.0      
## [33] rstudioapi_0.13    jquerylib_0.1.4    generics_0.1.2     jsonlite_1.8.0    
## [37] magrittr_2.0.2     Matrix_1.3-4       Rcpp_1.0.8         munsell_0.5.0     
## [41] fansi_1.0.3        lifecycle_1.0.1    stringi_1.7.6      yaml_2.3.5        
## [45] rootSolve_1.8.2.3  grid_4.1.2         crayon_1.5.0       lmom_2.8          
## [49] lattice_0.20-45    haven_2.5.0        hms_1.1.1          knitr_1.39        
## [53] pillar_1.7.0       boot_1.3-28        gld_2.6.4          codetools_0.2-18  
## [57] reprex_2.0.1       glue_1.6.2         evaluate_0.15      data.table_1.14.2 
## [61] modelr_0.1.8       vctrs_0.4.1        tzdb_0.3.0         cellranger_1.1.0  
## [65] gtable_0.3.0       assertthat_0.2.1   xfun_0.31          broom_0.8.0       
## [69] e1071_1.7-9        class_7.3-19       units_0.8-0        ellipsis_0.3.2